function vadata = genvadata(D,modspecs,controls)

if isstruct(D)
    if ~isequal(fieldnames(D),fieldnames(modspecs),fieldnames(controls))
        error('structure variables do not have same field name order')
    end
    vdat = cellfun(@(f) genvadata(D.(f),modspecs.(f),controls.(f)),fieldnames(D),'unif',0);
    teacherid = cellfun(@(x) reshape([x.teacherid],[],1),vdat,'unif',0);
    teacherid = unique(cat(1,teacherid{:}));
    vadata = struct('teacherid',num2cell(teacherid),'Y',[],'Z',[],'ZF',[],'sid',[]);
    vadata = repmat(vadata,[1 numel(vdat)]);
    o = 1;
    for f = fieldnames(D)'
        [vadata(:,o).Z] = deal(zeros([0 numel(modspecs.(f{1}).tchcoef)]));
        [~,idx] = ismember([vdat{o}.teacherid],[vadata(:,o).teacherid]);
        vadata(idx,o) = vdat{o};
        o = o + 1;
    end
    return
end

if ~any(strcmp(D.Properties.VariableNames,'intercept'))
    D.intercept(:) = 1;
end

D.row_id = (1:size(D,1))';
[uni_teacherid,~,D.tch_id] = unique(D.teacherid);
D = sortrows(D,{'tch_id','row_id'});

X = D{:,controls.Properties.RowNames};
Z = +D{:,modspecs.tchcoef};
Y = D{:,modspecs.outcome} - X*controls.Coefficient;
SID = D.sid;

[~,CHX] = fullcolrank(X);
if ~all(round(reshape(Z - CHX*((CHX'*CHX) \ (CHX'*Z)),[],1),8)==0)
    warning('Teacher coefficients not included in controls')
end

rep_tch = accumarray(D.tch_id,1);
if ~all(cellfun(@(j_id,j) all(j_id==j),mat2cell(D.tch_id,rep_tch,1),num2cell((1:numel(uni_teacherid))')))
    error('Problem splitting teacher data')
end

Y = mat2cell(Y,rep_tch,size(Y,2));
Z = mat2cell(Z,rep_tch,size(Z,2));
[~,ZF] = cellfun(@(z) fullcolrank(z),Z,'unif',0);
SID = mat2cell(SID,rep_tch,size(SID,2));

vadata = struct('teacherid',num2cell(uni_teacherid),'Y',Y,'Z',Z,'ZF',ZF,'sid',SID);

if ~all(cellfun(@(i) all(diff(i)>0),{vadata.sid}))
    error('sid data is out of order')
end

